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ABSTRACT 

We formulate a perturbative approximation to gravitational instability, based on La- 
grangian hydrodynamics in Newtonian cosmology. We take account of 'pressure' effect 
of fluid, which is kinematically caused by velocity dispersion, to aim hydrodynamical 
description beyond shell crossing. Master equations in the Lagrangian description are 
derived and solved perturbatively up to second order. Then, as an illustration, power 
spectra of density fluctuations are computed in a one-dimensional model from the 
Lagrangian approximations and Eulerian linear perturbation theory for comparison. 
We find that the results by the Lagrangian approximations are different from those 
by the Eulerian one in weakly non-linear regime at the scales smaller than the Jeans 
length. We also show the validity of the perturbative Lagrangian approximations by 
consulting difference between the first-order and the second-order approximations. 

Key words: gravitation - hydrodynamics - instabilities - cosmology: theory - large- 
scale structure of Universe. 



1 INTRODUCTION 

It is significant to investigate evolution of inhomogeneities by gravitational instability in the expanding universe from the 
viewpoint of cosmological structure formation. In order to find how to form cosmic structures via gravitational instability, 
numerical simulations such as A'^-body simulations have been carried out by several groups (e.g. Miyoshi & Kihara 1975; 
Hockney & Eastwood 1988; Couchman 1999). Such numerical approaches have brought us many useful informations about 
structure formation, but analytical approaches are also needed to obtain physical understanding of structure formation. 

For analytical approaches, one usually treats matter contained in the universe as a self-gravitating fluid, and considers 
solving the hydrodynamical equations for the fluid. Since the hydrodynamical equations are generally non-linear, one cannot 
solve them without any assumption or approximation. A conventional approximation is linear perturbation in homogeneous 



and isotropic universes, based on the Eulerian picture of hydrodynamics (Weinberg 1972; Peebles 1980; Coles & Lucchin 199. 



Sahni & Coles 1995). This approach is, by construction, valid only in linear regime, where amplitude of density perturbations 



is much smaller than unity. For description beyond linear regime, Zel'dovich (197C) proposed a new approximation scheme 



in which perturbations are given as the Lagrangian displacement of fluid flow by an extrapolation of the linear perturbation 
theory. This approximation is known as Zel'dovich approximation, which has been found to give relatively accurate results 



and work better than the Eulerian approximations by comparison with exact solutions (Munshi, Sahni & Starobinsky 1994 



Sahni & Shandarin 1996, Yoshisato, Matsubara & Morikawa 1998), in weakly non-linear regime, where amplitude of density 



perturbations becomes comparable to unity. The Zel'dovich approximation is shown to be a subclass of the flrst-order solutions 
of a perturbation theory in the Lagrangian hydrodynamics (Buchert 1989, 1992), and along this line, higher-order extensions 
have been developed up to third order (Bouchet et al. 1992; Buchert fc Ehlers 1993| ; Buchert 1994; Bouchet et al. 1995; [Catelan 



1995 



Sasaki fc Kasai 19981) 



In the Zel'dovich approximation, however, physical singularities called 'shell crossing' inevitably occur. This is a conse- 
quence due to the fact that a self-gravitating pressureless fluid is taken as a matter model in the approximation. At the epoch 
beyond shell crossing, the Zel'dovich approximation soon becomes inaccurate because the fluid elements move throughout in 
the directions which are set initially, and then inhomogeneous structures, which are formed compactly once, are dissolved in 
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the approximation scheme. To resolve this problem, some modifications have been proposed, such as 'truncated Zel'dovich 



approximation' (Coles, Melott & Shandarin 1993; Melott, Pellman & Shandarin 1994), which is an optimization of the ap- 



proximation by truncating small-scale fluctuations, and 'adhesion approximation' (Gurbatov, Saichev & Shandarin 1989), 
where an artiflcial viscosity is introduced. These modifications eliminate the shortcomings of the Zel'dovich approximation, 
and actually the modified approximations provide excellent results compared with A'^-body simulations in some cases. However 
the physical grounds of the modifications are not clarified. 

To have more well-founded approximations from a physical point of view, we need to study gravitational instability of 
pressureless matter beyond shell crossing. Buchert & Dommguez ( [L998D have examined this issue, starting from the collisionless 
Boltzmann equation, which is usually applied to the stellar systems (e.g. Binney & Tremaine 1987). They argued that effect of 
velocity dispersion will be significant beyond shell crossing, and if the velocity dispersion is approximately isotropic, it yields 
pressure-like or viscosity terms. This implies that the gravitational instability of pressureless matter beyond shell crossing 
can be described effectively by hydrodynamic equations for a fiuid with pressure-like force. Following this view, Adler & 
Buchert (1999) have proposed reformulation of the Lagrangian perturbation theory by taking account of pressure effect. They 
derived first-order perturbation equations in the Lagrangian coordinates under the assumption that the pressure is a function 
of only mass density. They did not, however, present solutions of the perturbation equations or analyse the evolution of 
density perturbations with the solutions. One may expect the reformulation to extend the regions which can be described 
by analytical approximations, but this should be confirmed by a concrete illustration. The aim of this paper is to show how 
the reformulation gives description of gravitational instability by solving perturbation equations and illustrating behaviour of 
density perturbations. 

In this paper, we derive and solve the Lagrangian perturbation equations with pressure up to second order, assuming a 
polytropic equation of state. We adopt the method of the Fourier transformation for the solutions and then will see mode 
couplings in the Lagrangian Fourier space in the second-order solutions. In particular, we obtain explicit form of the second- 
order solutions in the case P oc p*'^'^, where P and p are pressure and mass density, respectively. Moreover, as an illustration 
of the formulation, power spectra of density fiuctuations are computed in a one-dimensional model from the Lagrangian 
approximations for the case P cx p*''^ , and are compared with the results by the Eulerian linear perturbation theory to clarify 
the difference between them. We also compare the first-order and the second-order Lagrangian approximations to examine 
the validity of the approximation scheme. 

This paper is organized as follows. In Section 2 we present basic equations of our method. Starting from the hydrody- 
namical equations, we derive master equations of the Lagrangian perturbation theory with pressure effect. In Section 3 we 
obtain perturbation equations by expanding the master equations up to second order, and solve them via the Fourier trans- 
formation. Section 4 gives illustrative examples of computation of density perturbations by the Lagrangian and the Eulerian 
approximations in a one-dimensional model. Showing power spectra of density perturbations, we discuss differences among 
the approximations. Section 5 contains concluding remarks. 



2 BASIC EQUATIONS 

We begin with basic equations of cosmological hydrodynamics for a self-gravitating fluid with energy density p and 'pressure' 
P. In coordinates x = r/a comoving with cosmic expansion, they are 

| + 3^p+iv..(p«) = 0, (1) 
dv d 1 . „ . 1 „ ^ 

— + -v + -{v-\7:,)v^g V^P, 2) 

ot a a pa 

X g = , (3) 

■ g = -47rGa(p - pb) , (4) 

where a = a{t) is the cosmic scale factor, pb — ph{t) is energy density of a homogeneous and isotropic (background) universe, 
and V and g represent respectively velocity fleld and gravitational field strength due to presence of inhomogeneity, and 
then may be called as 'peculiar velocity field' and 'peculiar gravitational field.' The 'pressure' we take into account here is 
kinematical one due to occurrence of velocity dispersion beyond shell crossing of dust flow, as stated by Buchert & Dommguez 
(1998), rather than thermodynamical one. Thus equation (^ is close to the Jeans equation, which is gained by taking moments 



of the collisionless Boltzmann equation (e.g. Binney & Tremaine 1987). 

In Lagrangian description of hydrodynamics, using the time derivative along the fluid flow 

d d I, „ , 

t: = 7i7 + - ^ ■ V- . 

at ot a 

equations (^ and (^) are rewritten as 
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^ + 3^P+^(V.-«) = 0, (5) 
at a a 

dv a 1 „ ^ 

— + -v=g V,P. (6) 

at a pa 

The coordinates x of trajectories of tlie fluid elements are expressed by Lagrangian coordinates q, defined by initial values of 
the coordinates x, in the form 

x = q + s(q,t), (7) 

where q and s represent the background Hubble flow and deviation of the flow from the background, respectively. The 
continuity equation (P) is then exactly solved as 

P = pbJ^^ , (8) 
or equivalently for density contrast 5 = {p — ph)/pb, 

(9) 

where J = det{dxi/dqj) — det{Sij + dsi/dqj) is the Jacobian of the transformation x ^ q. The peculiar velocity is written 
by definition as t; = as , and from equation the peculiar gravitational field becomes 

\ a a'' dp 

where an overdot (') denotes d/dt. Note that the square of the 'sound speed,' dP/dp, is a function of p, and can be written 
in terms of s by using equation if an equation of state is provided. Now all physical quantities are found to be written in 
terms of s and it remains only to find solutions for s. We obtain the following master equations for s from equations and 

(|)^ 

2-s)=Q, (10) 
a J 

V, • f s + 2-s - -L^ J-V, J I = -AnGpbir^ - 1) . (11) 
\ a a-^ ap J 

The relation between and = d/dq is obtained from equation (^) as 

dqi dqi dxj y dqt J dxj dxi dqt dxj 
Using this equation iteratively, we have 
d _ d dsj d _ d dsj f d dsk 9 \ _ d dsj d dsj dsk d 



" J " " " ^ (13) 
dxi dqi dqi dxj dqi dqi \dqj dqj dxk J dqi dqi dqj dqi dqj dxk 

The treatment is fully non-linear and exact so far. Combining equations (|lo|), (pT|), and (|l^), we can obtain perturbative 
solutions for s up to any order in principle. It should be emphasized that density p is treated non-perturbatively because of 
equation (tel), even if solutions for s are obtained in a perturbative manner. 



3 PERTURBATIVE APPROACH IN LAGRANGIAN COORDINATES 
3.1 Derivation of perturbation equations 

Let us proceed a perturbative approach in the Lagrangian description. We write the displacement vector s in a perturbative 
form s = s^-^^ + s^'^^ + • • •. Superscripts (1) and (2) denote first-order and second-order quantities in perturbative expansion 
with respect to amplitude e of primordial fluctuations. We make the perturbative expansion only for s, and p is not expanded. 
Then we may expect relatively accurate description for p by this formulation, even in non-linear regime. In the perturbative 
expansion, equation (^o|) gives to first order, 

V, X (s«+2^s(^') =0, (14) 
and to second order, 

, X (s(^) + 2^s(^))] ^ = u,ks^^ + 2^iW ) , (15) 

where (•),i denotes d/dqi. 

Next we consider equation (llll). The Jacobian J is expanded as 
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J = l + Vq-s(''+V, -s 



= (2) 



(V, 



fl)x2 _ (1)11) 



and then the square of the 'sound speed,' dP/dp, can be written as 

j2 , 



dP , , dP, , d^P, , „ (I) 2s 

Thus we obtain to first order, 
f dP 



a a-' dp 



(1)- 



47rG'pbVq ■ s 



(1) 



(16) 



and to second order. 



dp 



1 d^P / (1) 2 (1 



4-KGph 



dp 

ifs(l))2_i (1) (1) 



(17) 



In order to solve the perturbation equations, it is convenient to decompose s'^' and s^^-* into the longitudinal and the transverse 
parts in the form 



= (1) 



= VnS + S' 



= (2) 



where S and ^ are respectively first-order and second-order scalar functions, and S'^ and (^'^ satisfy Vq ■ = 0, Vq • (^""^ = 0. 
To note their physical meanings, the first-order longitudinal and transverse parts are related to linear density and vortical 
perturbations, respectively. At the second-order level, however, such a simple interpretation of the perturbation modes does 
not hold any more. The first-order perturbation equations ( |l^ and ( |l6| ) then become 



Vq X ( 



a 



Vl ( S + 2-S-47rGpbS - -^^(pb)V^S ] =0. 
\ a a2 dp J 

Under some adequate boundary conditions, these can be reduced as 



+2-S' =0, 

a 

5 + 2-S - 47rGpbS - -!-^(pb)V^5' = , 
a dp 



(18) 
(19) 

(20) 
(21) 



which are obtained by Adler & Buchert ( 1999 ). 

The second-order perturbation equations ( |l5| ) and ( [l^ ) are also rewritten in terms of the longitudinal and the transverse 
parts. Equation (^) reads 

1 dP, 



Vq X (C + 2-0 



a2 dp 



{Ph)^ijk{S^£j + Si ,j)\7 gS^ke ■ 



(22) 



The curl of this equation gives 



a 



(23) 



where is a source term, which is quadratic with respect to the first-order perturbations, of the form 



Qi (q,i) 



47rGpb [Sj AkSjk + Sj ^iVqSj — VqSj S ,ij — Sj ,kS,ijk ) 



1 dP, 



+ ^-jJ^(Pb) {Sjjk'VqSjk + 'S',ij Vq VqS.j — VqS'jVqS.ij — SjfcVq5,;jfc) 



+ 



_l_dP 

a2 dp 



(Pb) (Sj ^ik\7qSjk + Sj ^i\/q\/qSj — VqSj VqS^ij — S j ^kV qS ^ijk') 



Equation (|17[) becomes 



Vq ( C + 2^C-4^GpbC-^^(/5b)V^C) =0(q,t): 



where 

Q{q,t) 



27rGpb [S,ijS,,j - [Vlsf] 



(24) 
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1 d P 1 P 

1 dp 

dp 



(Pb) (S^ijWlsJj + 2S^ijkSj^jk + vISa'^IsJ + 2VqS,ij5'7'j) 



—2nGpbSi'_jSj^i — ^"JJ^(/'t>) {sJ^jVqSjj + SjjkSjjk) ■ 

We can easily co nfirm that equations (^3|) and (^ ) are consistent with the second-order perturbation equations obtained by 
Sasaki & Kasai (1998) for the pressureless case. 



3.2 Solutions of perturbation equations 

Here we solve the perturbation equations in the presence of pressure effect. We assume that the background universe is 
spatially flat one with a{t) = t^^^ and pb = l/{6nGt^) oc for simplicity. The first-order perturbation equation ( po| ) for the 
transverse part has the same form as in the pressureless case. Thus we immediately find the solutions 

S'^ oc const. , t~'^^^ . (25) 
For the longitudinal part, the Fourier transform of equation ( |2l| ) with respect to the Lagrangian coordinates yields 

5 + 2-S-47rGpbS+-!-^(pb)lK|'5 = 0, (26) 
a dp 

where (•) denotes a Fourier component, and if is a wavenumber vector associated with the Lagrangian coordinates q. It 
should be noted that the form of equation (^) is similar to that of an equation for the density contrast Sun in the Eulerian 
linear theory. Actually it reads 

+ 2 7^-^-47rGpb<5iin fc,t -H— — Pb fc^5,i„ fe,t =0, 27 

ot^ a at dp 

where k denotes a wavenumber vector associated with the Eulerian coordinates a;. As an example, assuming a polytropic 
equation of state P — np"' with a constant k and a polytropic index 7, the solutions of equation ( p7[ ) are (Weinberg 1972) 

r t^^'^Ji, (A\k\t-''+^''') for 7 / I , 

^-(^'^)-{ ,-i/«±V5^7itW ^ for 7!}; 

where = 5/(8 — 67), J±v is the Bessel function of order ±v, and 

V^7(67rG)(^ R_4,„„_i/3 
^= ' B=-.{%.G) I . 

Note that the above solutions include wavenumbers of fiuctuations, whereas the solutions for the pressureless matter do not. 
We then find solutions of equation ( p^ with the help of the known results for the density contrast in the Eulerian linear 
theory. Hence, in the case of the polytropic equation of state P = Kp^ and if = 5/(8 — 67) is not an integer, we obtain 
general solutions for S{K, t) as 

S{K, t) =D+{K, t)C+ (K) +D-{K, t)C' (K) , (29) 
where D^{K,t) are provided, by replacing k with K in equation (|28|), in the form 

and C^(K) are determined by initial conditions. Notice that K is the Lagrangian wavenumber, which is different from the 
Eulerian wavenumber k. 

Next we consider solutions of the second-order perturbation equations ( ^ ) and (ji^). The Fourier transform of equa- 
tions ^) and (§^ gives 

\K\' (C + -Q^iK,t), (31) 



a 



1 dP, ,,_,2; 



- \K\' C + 2-C-4^GpbC + --^(Pb)!Ji:rC =QiK,t). (32) 
\ a a'' dp J 

The solutions are formally written as 

C{K,t) = -^[ dt'G'^{t,t')^{K,t'), (33) 
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\K\'^ 



dt' G{K,t,t')Q{K,t') 



(34) 



by using the Green functions G'^{t,t') and G{K ,t,t'). The Green function G'^{t,t') does not depend on an equation of state 
and is given as 



(35) 



while the G{K,t,t') depends on an equation of state. Under the assumption P = up'' , if 7 7^ 4/3 and v — 5/{8 — 67) is not 
an integer, we have 



G{K,t,t') = 



and if 7 = 4/3, 



(--I)"' 



^-1/6^/7/6 



2 sin u-K 

-J,iA\K\t-''+''/'')J-,{A\K\t'-''+-'^'-^ 
-1/2 



J^4A\K\t-'-'+^^'')J4A\K\t'- 



7+4/3-, 



(36) 



G{K t t') = ~- ^— - -B|X t'^^'^t''^^^ ^^-V25/36-S|K|2^/y25/36-Slif [2 _ ^y/2h/3a~B\K\i _i_l-^2i,/36-B\K\^\^ ^gr,-^ 

In order to present exphcit form of the second-order solutions, we must compute the Fourier-transformed source terms and 
Q. Hereafter we neglect the first-order transverse part S'^ in and Q for simplicity. This is equivalent to give no attention 
to effect of vorticity in the second-order solutions. They are written in the following convolution form; 



QT(K,t) 



i 1 AP , 



(pb) / A?K' S{K',t) S{K - K',t)\K - K'\^ K' ■ {K - K') 



(27r)3 a2 dp 

K' [K' ■ {K - K')) + K'\K - K'\^ - {K - K')\K'\'' - {K - K') [K' ■ {K - K')) 



(38) 



Q{K,t) = 



1 






^/ 


1 




+ — 


^< 


1 

+ — 


d^P 


dp2 



d^K' S{K',t) S{K-K',t) 



27rGpb 



(K' ■ {K - K')Y - \K'\^\K - K'\ 



\K'\'\K - K'\^ K' {K- K') + [K' ■ [K - K')Y + 2\K - K'\^ [K' ■ {K - K'))' 



(39) 



Using the first-order solution (|29|), we obtain 

C{K,t) = -^^J d'K'E''{K,K',t)ic+{K')C+{K^K')+C+{K')C-{K^K') 

+C~{K')C+{K - K') +C~{K')C'{K - K')^ \K -K'\^ K' [K - K') 
K' {K' ■ [K - K')) + K'\K - K'f - {K - K')\K'f - {K ~ K') {K' ■ {K - K')) 

^^^'^^ = ' (2^)3 d'K'(c+{K')C+{K - K') + C+{K')C-{K - K') 

+C-{K')C+{K - K') + C-{K')C-{K - K')^ E{K,K',t) [K' ■ {K - if'))' - I^^Tl-^^ - K' 
+Fi{K,K',t)\\K'f\K - K'f K' ■ {K - K') + {K' ■ {K - K'))'' + 2\K - K'f {K' ■ {K - K'))' 
+F2{K,K',t) [\K'\'^\K - K'f + \K'f\K - K'f K' ■ {K - K')] 



(40) 



(41) 



where time-dependent factors are given as 
a2(f') dp 

+D'{K',t')D+ {K - K',t') + {K',t')D' {K - K' , f 



E'^{K,K',t) = j Jj^^{pi,{t'))G'^{t,t') {D+(K',t')D+iK - K',t')+D+iK',t')D-{K - K',t') 



(42) 



© 2001 RAS, MNRAS 000, §-0 



Extending Lagmngian perturbation theory to a fluid with velocity dispersion 7 

E{K,K',t) = j dt'2nGpi,{t')G{K,t,t')(^D+{K',t')D+{K-K',t') + D+{K',t')D~{K-K',t') 

+D- {K' ,t')D+ {K - K' ,t') + D- {K' ,t')D- {K - K' ,t')^ , (43) 

F,{K,K',t) = J -^^{p^{t'))G{K,t,t'){D+{K',t')D+{K~K',t')+D+{K',t')D-{K~K',t') 

+D- {K',t')D+{K - K',t') + D- [K', t')D- [K - K',t')j , (44) 

F2iK,K',t) = / -^^(pUO)pb(^0G(lf,^,^0(o+(-^f^Oo+(-^^--^f^^0 + ^+(J^^^0^"(J^-^^^O 

+D' {K',t')D+{K - K',t') + D' {K', t')D- {K - K' ,t') \ . (45) 



The convolution in the solutions ^ and (0) represents mode couplings in If-space, which inevitably occur at second order 
due to non-linearity. Although it is cumbersome to perform the integration in equations (^)-(^^ in general, it is easy to do it 
if an equation of state is P = np^^^ . In this case, if only the D'^ {K' ,t)D'^ {K — K' ,t) part is considered, equations (|42[)-(|4^) 
become 

f,,-l/3+b(K')+b{K~K') 

E^{K,K',t) = — , (46) 

[^\ + h{K')+h{K-K')\[h{K')+h{K-K')] 

,-l/3+b{K')+b(K-K') 

E(K,K',t) = ^, (47) 

3[-^-b{K) + biK') + b{K^K')][-l + b{K)+b{K')+b{K-K')] 

I j^-^-l/A+b{K') + b{K-K') 

F^[K,K ,t)= [^l+h{K) + b{K') + b{K-K')\ ' ^^^^ 

r,f--l/'i+b{K')+b{K-K') 

Fo{K K' t) = — (49) 

i[-l~b{K)+b{K')+b{K-K')\[-l + b{K) + b{K')+b{K-K')\' 

where b{K) = ^25/36 - B\K\^. Note that all these factors have the same temporal dependence, f-i/3+Kif ')+i>{if -Jf ') Qf 
course, it is not a general property of the second-order solutions, but F\{K ,K' ,t) and F2{K , K' ,t) always have the same 
temporal dependence as long as an equation of state is of the form P — Kp^ . 



4 ILLUSTRATION IN A ONE-DIMENSIONAL MODEL 

In this section, we present illustrative examples of computation by the Lagrangian perturbation theory formulated in the 
previous section. We compute power spectra of density perturbations by the Eulerian linear theory, the Lagrangian first- 
order and second-order approximations in a one-dimensional model and then clarify difference between the Eulrian and the 
Lagrangian approximations. In linear regime \5\ <C 1, the Eulerian and the Lagrangian approximations give the same results. 
However, when these approximations are extrapolated into non-linear regime, the results given by them do not always coincide. 

In the pressureless case, the first-order approximation (i.e. the Zel'dovich approximation) coincides with an exact solution 
in a one-dimensional model. Although this does not hold in the presence of pressure in general, we discuss that the Lagrangian 
perturbative approximations may provide nearly exact description in weakly non-linear regime \5\ ^ 1, by consulting difference 
between the first-order and the second-order approximations. Of course, we cannot say that one-dimensional examples are 
realistic, but they are instructive to show advantages and features of non-linearity which the Lagrangian perturbation theory 
involves. 

4.1 Equations and perturbative solutions in a one-dimensional model 

First we present basic equations and perturbative solutions in a one-dimensional model. In the Eulerian linear approximation, 
the density contrast satisfies equation (^^, that is 

'-"^ ' + 2 ""\ - 47rGpb(5ii„ fc, t) + —— pb fc fc, t = , 50 

ot^ a at dp 

where k is the first component of the Eulerian wavenumber vector k. We find from this equation that density perturbations 
whose wavenumbers are smaller than 
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dP/dp 

will grow to form inhomogeneous structures, and those whose wavenumbers are larger than fcj will decay with acoustic 
oscillations (the Jeans condition). In particular, we immediately see the behaviour of density perturbations in the case 
P — Kp^l^ , where the solutions of equation are 



5un(fc,t)0Ct-^/«±V25/36-Bfe2^ ^^^^ 



and fej = ^2/(3B). If fc < fcj, one of the solutions becomes a growing solution, whereas if A; > fcj, both of the solutions are 
decaying ones. 

The relation between the Eulerian and the Lagrangian coordinates, equation (Q), in a one-dimensional model can be 
written as 

x\ = gi + si(iji,t) , 

X2 = 92 , (52) 

and energy density, equation (|^, is then 

p{qut)= J"'^'] (53) 
1 + si,i(gi,t) 

because the Jacobian J = 1 + S14. The relation between Va, and Vq, equation ([l^, becomes 

did d d d d , , 

dxi J dqi ' dx2 dq2 ' 8x3 dq-j, 



Hence from equation (|11|) we have 

(s\ + 2-s\ - 4^ J,i) = -4^Gpb(J-i - 1) , (55) 
J oqi y a a'' dp J 

which can be reduced as 

s'l + 2-si - 47rG'/0bSi - ^^'^^ =0. (56) 
a dp (1 + Si,i)^ 

by imposing appropriate boundary conditions. If we assume that an equation of state is of the form P = Kp"*", we obtain by 
using equation (|53|), 

,^i+2^si-4^Gpb.i-^^^— ^iii-^=0. (57) 
a (l + si,i)i+T 

It seems to be difficult to solve equation (ji^) or equation ( ^tI) exactly in general, although Gotz ( 1988 ) solved it in the 
case P (X p without cosmic expansion. Then we consider their solutions in a perturbative manner and adopt the perturbation 
solutions obtained in the previous section. Note that the displacement vector s — {si, 0, 0) consists only of the longitudinal parts 
{S,i, • • ■), because the transverse parts {S'^ , ^'^ , . . .) vanish in a one-dimensional model. The perturbation solutions ( |29| ) 
and (|lt) become 

SiK,t) = D+{K,t)C+{K) + D'{K,t)C'{K) , (58) 



(iK,t) = -TT^ / dK' (C+{K')C+{K-K') + C+{K')C-{K-K') + C-{K')C+{K~K') + C-iK')C~{K~K')) 
ItvK f ^ ' 

■ [2Fi{K,K',t) + F2{K,K',t)]K'^(K - K'f , (59) 

where K is the first component of the Lagrangian wavenumber vector K. The part proportional to E(K , K' ,t) in the 
second-order solutions ( ^H ) does not appear in the above expression because it vanishes in a one-dimensional model. 

In order to simplify the perturbation solutions further, let us consider the case where an equation of state is P = Kp*/^. 
Although the validity of this assumption is not clarified, it would be useful to understand features of the perturbation theory 
we have formulated. The temporal factors are computed as 

D^{K,t)^t-^^''^''^''\ (60) 

7 T,.-l/3+b(K') + b(K-K') 

2Fi(K,K\t) + F2(K,K',t) ^ --^^ ^, (61) 

S[^^-b{K)+b{K')+b{K^K')][^^ + b{K) + b{K') + biK-K')] 

where we again take into account only the D'^ {K')D^ [K — K') part to compute 2Fi + F2. 
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4.2 Initial conditions 



We consider setting of initial conditions for illustration. Initial conditions for two independent physical quantities are needed 
to determine C^{K). Here we impose initial conditions for the density contrast S and the peculiar velocity v, whose initial 
values are denoted by Sin and Uin, respectively. For comparison with the pressureless case, we take Sin and vin as those given by 
the Zel'dovich approximation, which is a subclass of the first-order approximation for pressureless fluid and becomes an exact 
solution in a one-dimensional model. By this setting, the C^{K) are expressed in terms of only Sin, because the Zel'dovich 
approximation includes just one arbitrary spatial function, through which we can make a relation between Sin and Din. The 
one-dimensional Zel'dovich approximation is written as 

X2 = 92 , (62) 
X3 = 93 , 



P(9i,i) 



where ^'(91) is an arbitrary spatial function, describing initial inhomogeneity. From these equations, we have 



Sin{qi) = 



l-ft2/3^_ll(qi) 



Viniqi) = ((2/3)t^''^*,i(gi),0,0) 



11(91), 



= ((2/3)*,i(gi),0,0) 



1. On the other hand, the first-order solution in the case P - 



where we define an initial time tin 

K{K)=K\C+{K) + C-{K)), 

i^,{K) = {iK [{-1/6 + b{K))C+{K) + [-1/6 -b{K))C-{K)] ,0,0) . 
Comparing equations (p^, (p5[), (p6[), and (p7[), we obtain 
5^{K) [ ^ . 5 



Kp*^^ gives 



2A'2 



1 ± 



6fe(A') 



(63) 

(64) 
(65) 

(66) 
(67) 

(68) 



Thus C (K) are completely determined if 5in{K) is provided in some appropriate manner. In our illustration, we choose 
Sin(K) = \Sin{K)\ exp(i(j)K) so that |(5in(-?s')|^ oc l-K^I", where n is a spectral index, and the phases (f>K are randomly distributed 
on the interval [0,27r]. This choice is a simpliflcation of the Gaussian statistics for initial density perturbations Sin{qi) in real 
space (Coles & Lucchin 1995), which is usually adopted in the study of large-scale structure formation. 



4.3 Evolution of power spectra of density perturbations 

Now we compute power spectra P(|A;|, t) = P(k, t) = {\S{k, f)p) of density perturbations, where (•) denotes ensemble average 
over the entire distribution, by using the Lagrangian perturbation theory formulated in the preceding section. We also 
compute them by the Eulerian linear theory and compare the results to clarify difference of the Eulerian and the Lagrangian 
perturbative approximations. In the Lagrangian approximations, we need some computation to obtain the power spectra 
although the Eulerian approximation yields them directly. The procedure of the computation is the following: 

(i) First we specify initial conditions as we mentioned above. Then we have complete form of the perturbation solutions in 
if-space. 

(ii) Next we transform the perturbation solutions in Tf-space into those in g-space via the inverse Fourier transformation. 

(iii) From the perturbation solutions in g-space, we immediately flnd density perturbations in g-space by equation (^3|). 

(iv) We evaluate density perturbations in a;-space from those in g-space. 

(v) Finally by the Fourier transformation with respect to x, we obtain power spectra of density perturbations. 

By way of this procedure, we obtain the power spectra of density perturbations presented in Figs |^-^ We choose a spectral 
index as n = -1-1, 0, and —1. The constant B, which is proportional to k and then provides strength of pressure effect, is put 
by hand. Here it is chosen so that the Jeans wavenumber fcj = y^2/{3B) is 80. Note that is now a constant because of 
our choice of the polytropic index 7 = 4/3, although fc,i depends on time in general. We set initial conditions at a = 1, and 
pursue the evolution up to a = 1100. 

Indeed, for all the cases n = -1-1, 0, and —1, shell crossing will occur at a ~ 1100 in the Lagrangian approximations despite 
the presence of pressure effect. (In other words, we normalize amplitude of initial density fluctuations so that shell crossing 
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- Eulerian linear 

- Lagrangian n=l 

- Lagrangian n=0 
-Lagrangian n=-l 




0.01 r 



0.0001 



Figure 1. The 'transfer function' of density perturbations at a = 300 
computed from the Eulerian linear theory and the Lagrangian first- 
order approximations. It does not depend on the initial conditions in 
the Eulerian linear theory, but does in the Lagrangian approxima- 
tion. 



Eulerian linear 

— Lagrangian n= 1 

Lagrangian n=0 

Lagrangian n=-l 




0.01 r 



0.0001 



Figure 2. Same as in Fig. ^, but for a = 1100. 




Figure 3. Power spectra of density perturbations at a = 1100 
computed from the Eulerian linear theory, and the Lagrangian 
first-order and second-order approximations. The spectral in- 
dex is n = — 1. 



Figure 4. Power spectra at a = 1100 for n = —1, showing 
small difference between the results by the first-order and the 
second-order Lagrangian approximations. The spectrum by 
the Eulerian linear approximation is omitted. 



occurs at a ~ 1100 in the Lagrangian approximations.) Then we cannot follow the evolution so deeply into non-linear regime, 
as long as we consider the stages before the occurrence of shell crossing. Actually in the illustration, "Pd^D ~ 10~'*-10~^ at 
a — 300, and V{\k\) ~ 10~"^-10~'' at a = 1100. Hence we must say that density perturbations remain in weakly non-linear 
regime (or in linear regime), rather than in non-linear regime, through this illustration. 
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Figs |l| and ^ show our results obtained by the Eulerian hnear theory and the Lagrangian first-order approximation, in 
terms of the 'transfer function,' 'P{\k\,t)/'Pin{\k\). It is convenient to use it because it does not depend on the initial conditions 
in the Eulerian linear theory, which actually yields 

2 

(69) 



Vlin{\k\,t) _ 1 



VU\k\) 4 

In the Lagrangian approximations, however, it generally depends on the initial conditions. We do not present in Figs |l| and ^ 
the results by the Lagrangian second-order approximation, because they are almost coincident with those by the first-order 
one. To show the difference between the Lagrangian first-order and second-order approximations, we show in Figs I and I the 
power spectra at a = 1100 for n = — 1, where the difference is largest within our calculations. 

First we observe the results by the Eulerian linear theory. The power spectra presented in Figs ^ and |^ show the very 
behaviour of linear density perturbations stated in subsection 4.1, i.e. density perturbations with wavenumbers \k\ < kj grow 
while those with wavenumbers |fe| > fcj decay with acoustic oscillation. On the other hand, in the Lagrangian approximations, 
the shape of the curves is manifestly different from that in the Eulerian linear theory. Although there is little difference on 
large scales (|fc| < fcj), we see that amplitude on small scales (|fc| > fcj) in the Lagrangian approximations is larger than that 
in the Eulerian linear theory. This fact has been observed also in the pressureless case ( Schneider & Bartelmann 1995 ) . The 
difference is small at a = 300, but becomes larger as time proceeds. Comparing the Lagrangian first-order and second-order 
approximations in Fig. ^ they are found to give almost coincident results through our computation from a = 1 to a = 1100 
and we can hardly observe difference between them. We present the enlarged power spectra in Fig. ^ where the difference is 
barely visible. Indeed the difference is less than 10% at jfcj ^ 150. In the second-order approximation, however, amplitude of 
the power spectrum is slightly suppressed, compared with the first-order one. The features of the power spectra mentioned 
above are common for all the cases, n = -1-1, 0, and —1. 



4.4 Discussions on the power spectra 

Let us consider the reasons of the features of the power spectra. First we examine the difference between the Eulerian and 
the Lagrangian approximations, which has been discussed in the pressureless case by Schneider & Bartelmann ( |199E| ). In 
the Eulerian linear theory, density perturbations with a wave mode evolve without being influenced by those with another 
mode. In the Lagrangian approximations, however, non-linearity are induced in calculation of density perturbations. Origin of 
non-linearity exists in an expression of the density contrast in the Lagrangian description as well as in the transformation of 
the density contrast in the Lagrangian coordinates q into those in the Eulerian coordinates x. To see this fact, let us express 
the Lagrangian density perturbations in the Eulerian coordinates in a one-dimensional system within a perturbative manner. 
The relation between x and q is given as 

xi = qi + si{qi,t) , 

where we assume that si is small enough to be treated as a perturbation. Using this relation iteratively, the inverse relation 
is obtained as 

qi = xi — si{qi,t) = xi — si{xi,t) + O ((si)^) . 
Then we have 

si{qr,t)^sr{xi,t)~^^l^^s^{xi,t) + 0{{sif) , (70) 

d dxi d d dsi{qi,t) d d dsi{xi,t) d , , 

dq\ dqi dx\ dxi dqi dxi dxi dxi dxi 

By using equations (^o|) and ([ti]), the expression of the density contrast in the Lagrangian coordinates 

5,(g„t) = (1 + s,A,i,t)r ~ 1 = -^^^ + O {{s.f) 

(a subscript 'L' denotes 'Lagrangian') is transformed as 

dsi[xi,t) d^si{xi,t) ( dsi{xi,t) \^ , 

^^("^'') = dxi +0((^i)), (72) 

where the first term of the right-hand side corresponds to the density contrast 5ii-n{xi,i) in the Eulerian linear theory. 
Equation ([t^) indicates that the density contrast obtained by the Lagrangian description includes extra non-linear terms, 
which cause mode couplings. For example, if initial density perturbations consist of a single wave mode, 5i-n(x\) oc sinfca::i, such 
non-linear terms generate high-frequency modes such as sin 2kx\ . Of course, in the Eulerian linear theory, mode couplings never 
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occur and existence of just a single mode is preserved, i.e. 5iin(a;i,t) oc sinfca::i. To see the mode-coupling effect quantitatively, 
one should consider the Fourier transform of equation (^), 

S^ik,t) = 5,Uk,t) + ^ >^SUk',t)SUk k',t) ^ ^ ^^^^^^^3) (,3) 

where the second term of the right side represents the mode couplings. Note that the second term is written as the summation 
with respect to all wavenumbers. Then the effect of the second term may be larger than the simple square of Siin{k,t) by 
order of the number of wave modes. In the case of our illustration, it may be larger by 2 order. For example, if the power 
spectrum has a peak value 10~^, the mode-coupling effect can generate the amplitude of 10~® at high frequency. One may 
wonder, at first glance of our results, why the results by the Lagrangian approximations contain so large amplitude at high 
frequency, although the illustration is performed in nearly linear regime. However, this fact is explained by the effect of the 
mode couplings. 

The appearance of the large amplitude on small scales may be interpreted physically as follows. In the Lagrangian 
description of hydrodynamics, one obtains physical quantities in a frame comoving with flow lines of fluid. If there exists 
growing density enhancement in a region, one can see that flow lines there become close to each other because of gravitational 
instability. In other words, one knows by following flow lines that a physical wavelength of inhomogeneity gets small due to 
gravitational contraction. Actually, density perturbations with an initially small wavenumber \K\ become those with a large 
|fc| later. It can easily seen by the relation between K and k, given as 

' 1 A^MMV^i,,i(,+,(,,,)), 



' ' I L + si(L,t) L\ L 

where ^ is a physical wavelength of inhomogeneity measured in the Eulerian coordinates, and L denotes an initial wavelength. 
Thus we may conclude that the appearance of the large amplitude on small scales is due to the fact that scale of inhomogeneity 
is shortened as inhomogeneity grows because of gravitational instability. 

Furthermore, let us make sure of behaviour of the power spectra by the Lagrangian approximations for large |A;|. As we 
stated in the previous subsection, shell-crossing singularities arise in spite of the presence of pressure effect in our perturbation 
scheme. In our illustration, they arise at a ~ 1100, and thus the epoch a = 1100 is just before the occurrence of shell crossing. 
In such an epoch, the power spectrum behaves like Pd^l) oc \k\~^ for large |fc| in a one- dimensional system ( ^chneider fc 
Barte mann 1995). This behaviour concerns only the occurrence of shell crossing, and is seen not only in the Zel'dovich 
approximation but also in our results. Figs ^ and ^. Note that in Fig. ^ the 'transfer function' behaves like 'P(|fc|)/'Pin(|fc|) ~ 
at high frequency, showing dependence on initial conditions, while Fig. |^ shows the behaviour directly. It should be 
also stressed that the above three kinds of the arguments on the Lagrangian power spectra hold true, whether the pressure 



effect is taken into account or not. In this sense, it is natural that the results by Schneider & Bartelmann (1995) and ours 
have similar tendency. 

Next let us confirm little difference between the first-order and the second-order Lagrangian approximations. For a rough 
estimation, we consider perturbations with a single wave mode K so that the first-order solution in the g-space is written in 
the form 

Siqi,t) = -^D(K,t)smKqi, (74) 
where e denotes amplitude of initial density perturbations, and D{K, t) — ^-i/^+^'C^) Then the second-order solution becomes 
a<l^,t)^~^^^D{K,tfsin2Kq,, (75) 

where Kj = fc,i — a/2/(3B) is a wavenumber corresponding to the Jeans length. The fraction of S{qi,t) and ({qi,t) is 
estimated as 



C(9i,i) 






S{qi,t) 







K,t). (76) 

\i K > Kj, the factor D{K, i) decreases as time and then the fraction \C,/ S\ remains small forever. In contrast, li K <^ Kj, the 
factor D{K,t) increases as time, but K/Kj is small. Then the fraction \C^/S\ cannot grow to be so large in early time. Indeed 
in our calculations, it is less than about 10% during a period up to a = 1100. As time proceeds, however, it will become large 
ii K <^ Kj. This is a simple argument, but inequality (^^ may be useful to give a criterion of effect of the second-order 
terms. To let this argument more rigorous, we should take into account the mode-coupling effect, as we do in equation ([tsI). 
In the estimation of the second-order solution, however, it is not essential to include the mode-coupling effect. It is rather 
significant to notice that the right side of inequality (|7^) has the factor (K/Kj)^. The presence of this factor is due to the 
fact that the second-order solution is of purely pressure origin in a one-dimensional model. In other words, gravitational effect 
is completely included in the first-order solution in a one-dimensional model. Thus we can confirm small difference between 
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the first-order and the second-order Lagrangian approximations. This fact teils us that difference between the first-order 
approximation and an exact soiution is aiso small, at least up to a = 1100. On the other hand, in the pressureless case, the 
first-order approximation (i.e. the Zel'dovich approximation) becomes an exact solution in a one-dimensional model, and this 
fact is a strong ground of the validity of the Zel'dovich approximation. In this sense, we can also expect accurate description 
by the Lagrangian approximations in weakly non-linear regime in the presence of pressure, as in the pressureless case. 



5 CONCLUDING REMARKS 

We have developed a perturbative approximation theory, based on the Lagrangian description of hydrodynamics in the 



framework of the Newtonian cosmology, by extending the method of Adler & Buchert (1999). Including 'pressure' effect of 
fluid, we have derived and solved perturbation equations in the Lagrangian coordinates up to second order. Especially we 
presented explicit form of the second-order solutions for the case P oc p*^^. We have also computed the evolution of the power 
spectra of density perturbations in a one-dimensional model, based on the Eulerian and the Lagrangian approximations. 
Comparing the power spectra, we have found difference of these approximations. In particular, large amplitude on small 
scales has appeared in the results of the Lagrangian approximations beyond linear regime. Moreover, the first-order and the 
second-order Lagrangian approximations have been found to yield almost the same results within our calculations. Then we 
can conclude that, in a one-dimensional system, the first-order Lagrangian approximation provides nearly exact description 
in weakly non-linear regime. 

In the computation of the power spectra of density perturbations by the Lagrangian approximations, we have found 
that the shell-crossing singularities occur even in the presence of the pressure effect. However, the epoch of the occurrence of 
shell crossing in our approximations is, of course, later than that in the Zel'dovich approximation. This fact is easily seen by 
considering perturbations with a single wave mode so that the first-order solution is given by equation (Q) again. Then the 
energy density, equation (jHsj), becomes 

^fa'^)= l-.Z>(g,^tini..r 

If if < K}, the denominator of the right side goes to zero in a finite time, i.e. shell crossing will occur. But, since the growth 
rate D{K,t) is weaker than that of the Zel'dovich approximation, the shell-crossing epoch becomes later. 

In this paper, we focus our attention on the case P oc p*^^, where the solutions of the perturbation equations are written 
in a simple form. This equation of state is nothing but an assumption to simplify the perturbation solutions. However, it 
is crucial what form an effective equation of state takes when velocity dispersion is replaced with pressure-like force. Thus, 
in order to let our formulation more useful, we must reconsider an equation of state which holds effectively in high-density 



regions, where velocity dispersion plays an important role. For example, Buchert & Domfnguez (1998) found that a relation 
P oc p^^^ is favoured for small velocity dispersion under the kinematical restriction that the fiuid motion involves no shear. 
Extensions of our formulation to such cases, as well as more generic determination of an equation of state, will be the subjects 
of future investigation. 



ACKNOWLEDGMENTS 

We would like to thank the referee. Professor Bernard Jones, for constructive comments. We also thank Thomas Buchert, 
Kei-ichi Maeda, Hiroki Anzai, and Momoko Suda for helpful discussion and many valuable remarks. 



REFERENCES 

Adler S., Buchert T., 1999, A&A, 343, 317 

Binney J., Tremaine S., 1987, Galactic Dynamics. Princeton University Press, Princeton 

Bouchet F. R., Juszkiewicz R., Colombi S., Pellat R., 1992, ApJ, 394, L5 

Bouchet F. R., Colombi S., Hivon E., Juszkiewicz, R., 1995, A&A, 296, 575 

Buchert T., 1989, A&A, 223, 9 

Buchert T., 1992, MNRAS, 254, 729 

Buchert T., 1994, MNRAS, 267, 811 

Buchert T., Dommguez A., 1998, A&A, 335, 395 

Buchert T., Ehlers J., 1993, MNRAS, 264, 375 

Catelan P., 1995, MNRAS, 276, 115 

Coles P., Lucchin F., 1995, Cosmology: The Origin and Evolution of Cosmic Structure. John Wiley & Sons, Chichester 
Coles P., Melott A. L., Shandarin S. F., 1993, MNRAS, 260, 765 

Couchman H. M. P., 1999, in Miyama S. M., Tomisaka K., Hanawa T., eds. Numerical Astrophysics. Kluwer Academic Publishers, 
Dordrecht, p. 1 



© 2001 RAS, MNRAS 000, §-0 



14 M. Morita and T. Tatekawa 



Gotz G., 1988, Class. Quantum Grav., 5, 743 

Gurbatov S. N., Saichev A. I., Shandarin S. F., 1989, MNRAS, 236, 385 

Hockney R. W., Eastwood J. W., 1988, Computer Simulation Using Particles. lOP Publishing, Bristol 
Melott A. L., Pellman T. F., Shandarin S. F., 1994, MNRAS, 269, 626 
Miyoshi K., Kihara T., 1975, Pub. Astron. Soc. Japan, 27, 333 
Munshi D., Sahni V., Starobinsky A. A., 1994, ApJ, 436, 517 

Peebles P. J. E., 1980, The Large-Scale Structure of the Universe. Princeton University Press, Princeton 

Sahni V., Coles P., 1995, Phys. Rep., 262, 1 

Sahni V., Shandarin S., 1996, MNRAS, 282, 641 

Sasaki M., Kasai M., 1998, Prog. Theor. Phys., 99, 585 

Schneider P., Bartelmann M., 1995, MNRAS, 273, 475 

Yoshisato A., Matsubara T., Morikawa M., 1998, ApJ, 498, 48 

Weinberg S., 1972, Gravitation and Cosmology. John Wiley & Sons, New York 

Zel'dovich Ya. B., 1970, A&A, 5, 84 



© 2001 RAS, MNRAS 000, 



